In [1]:
# Definimos el cuerpo donde vamos a trabajar
m = 4;rango = 3;N = 2^m
K_.<a> = GF(2)
F.<a> = GF(2^m) #creado el campo donde "vivirá" el código
In [2]:
# Creamos el anillo de polinomios
PR = PolynomialRing(F,'X')
X = PR.gen()
g = X^3+X+1 # Polinomio de Goppa
L = [a^i for i in range(N)] # Soporte del codigo
In [3]:
# -------------------
# Funcion auxiliar que permite descomponer un polinomio en irreducibles
# -------------------
def descomponer_polinomio(p):
# El siguiente metodo permite descomponer un polinomio p en factores irreducibles p(z) = p0 (z) + z p1 (z)
# Entrada: Polinomio p
Phi1 = p.parent()
p0 = Phi1([sqrt(c) for c in p.list()[0::2]])
p1 = Phi1([sqrt(c) for c in p.list()[1::2]])
return (p0,p1)
# -------------------
# Algoritmo de euclides extendido: Obtener MCD y los s,t que lo generan.
# -------------------
def algoritmo_euclides_extendido(self, other):
delta = self.degree() #grado de polinomio 1
if other.is_zero(): # si el polinomio introducido es
ring = self.parent() #comprobamos el cuerpo en el que trabajamos
return self, R.one(), R.zero() #mcd = mismo polinomio y devuelve un uno (s) y un cero (t) en el cuerpo que trabajamos.
# mcd (a,b) = as+bt
ring = self.parent() #comprobamos el cuerpo en el que trabajamos
a = self # guardamos una copia del primer polinomio 1 (self)
b = other # guardamos una copia del segundo polinomio (other)
s = ring.one() # guardamos en s el uno del anillo
t = ring.zero() # guardamos en t el cero del anillo
resto0 = a
resto1 = b
while true:
cociente,resto_auxiliar = resto0.quo_rem(resto1) # La funcion quo_rem de Sage devuelve el cociente y el resto. Que guardamos en Q y ring.
resto0 = resto1
resto1 = resto_auxiliar
s = t
t = s - t*cociente
if resto1.degree() <= floor((delta-1)/2) and resto0.degree() <= floor((delta)/2):
break
V = (resto0-a*s)//b
coeficiente_lider = resto0.leading_coefficient() # guardamos el coeficiente lider del resto 0
return resto0/coeficiente_lider, s/coeficiente_lider, V/coeficiente_lider
# -------------------
# Funcion que calcula la inversa de un polinomio utilizando el algoritmo de euclides de Sage
# -------------------
def inversa_g(p,g):
(d,u,v) = xgcd(p,g)
return u.mod(g)
# -------------------
# Funcion de decodificacion de Patterson
# -------------------
def decodePatterson(y):
# Calculamos primero el vector alpha con los elementos primitivos.
alpha = vector(H*y)
# Consideramos nuestras matrices T,Y,Z definidas así como nuestro polinomio irreducible g
polinomioS = PR(0) # Inicializa como el polinomio 0 del anillo
for i in range(len(alpha)):
polinomioS = polinomioS + alpha[i]*(X^(len(alpha)-i-1)) # Lo vamos rellenando con los alpha
vector_g = descomponer_polinomio(g) # Guardamos en vector_g el par de polinomios irreducibles
w = ((vector_g[0])*inversa_g(vector_g[1],g)).mod(g)
vector_t = descomponer_polinomio(inversa_g(polinomioS,g) + X)
R = (vector_t[0]+(w)*(vector_t[1])).mod(g)
(a11,b11,c11) = algoritmo_euclides_extendido(g,R)
# Definimos el polinomio sigma
sigma = a11^2+X*(c11^2)
# Vamos comprobando uno a uno los coeficientes de sigma
# para asi determinar el conjunto de posiciones de error E - {i tal que e_i es distinto de 0}
for i in range(N):
if (sigma(a^i)==0): # an error occured
print ("Error encontrado en la posición: " + str(i))
y[i] = y[i] + 1
return y
In [4]:
T = matrix(F,rango,rango)
for i in range(rango):
count = rango - i
for j in range(rango):
if i > j:
T[i,j]=g.list()[count]
count = count + 1
if i < j:
T[i,j] = 0
if i == j:
T[i,j] = 1
In [5]:
print ("Matriz T: ")
show(T)
Out[5]:
In [6]:
V = matrix([[L[j]^i for j in range(N)] for i in range(rango)])
D = diagonal_matrix([1/g(L[i]) for i in range(N)])
H = T*V*D
H_Goppa_K = matrix(K_, m*H.nrows(),H.ncols())
for i in range(H.nrows()):
for j in range(H.ncols()):
be = bin(eval(H[i,j]._int_repr()))[2:];
be = '0'*(m-len(be))+be; be = list(be);
H_Goppa_K[m*i:m*(i+1),j]=vector(map(int,be));
show(H_Goppa_K)
Out[6]:
In [7]:
Krnl = H_Goppa_K.right_kernel();
G = Krnl.basis_matrix();
S = random_matrix(GF(2), N-m*rango)
while (S.determinant()==0):
S = random_matrix(GF(2), N-m*rango)
rng = range(N)
P = matrix(GF(2),N);
for i in range(N):
p = floor(len(rng)*random());
P[i,rng[p]]=1; rng=rng[:p]+rng[p+1:];
G_prima = S*G*P
show(G_prima)
Out[7]:
In [8]:
u = vector(K_,[randint(0,1) for _ in range(G_prima.nrows())])
c = u*G_prima
print 'Vector u' ; show(u)
print 'Vector c' ; show(c)
e = vector(K_,N)
e[8] = 1
e[9] = 1
print 'Vector de errores e' ; show(e)
y = c + e
print "Vector codificado y" ; show(y)
Out[8]:
Out[8]:
Out[8]:
Out[8]:
In [9]:
yP = y*(P.inverse())
yd = decodePatterson(yP)
corregido = (G.transpose()\yd)*S.inverse()
show(corregido)
Out[9]: